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A highly efficient Monte Carlo method for the calcula- 
tion of the density of states of classical spin systems is pre- 
sented. As an application, we investigate the density of states 
Qn(E, M) of two- and three-dimensional Ising models with N 
spins as a function of energy E and magnetization M. For 
a fixed energy lower than a critical value E Ci n the density of 
states exhibits two sharp maxima at M = ±M 3p (E) which de- 
fine the microcanonical spontaneous magnetization. An anal- 
ysis of the form M sp (E) oc (E c ,oo — E) 13 ^ yields very good 
results for the critical exponent /3 £ , thus demonstrating that 
critical exponents can be determined by analysing directly the 
density of states of finite systems. 

PACS numbers: 64.60. Fr, 05.10.-a, 05.50.+q, 64.60.-i 



I. INTRODUCTION 

Substantial progress has been achieved in the under- 
standing of discontinuous phase transitions since thfi 
problem has been formulated in microcanonical termsEm. 
Till now the microcanonical analysis of continuous phase 
transitions has only been partly successful. On the one 
hand it is a great accomplishment that the typical fea- 
tures of symmetry breaking, i.e. the abrupt onset of the 
order parameter as the critical point is approached from 
above and the diverging susceptibilities, turn up already 
for finite systemso. This is in contrast to the canon- 
ical ensemble where singularities appear exclusively in 
the thermodynamic limitp. On the other hand it is an 
intriguing and disappointing fact that mean field values 
have been founcLfor all critical exponents and for all fi- 
nite system sizesa when directly analysing the density of 
states. 

We do not contest the earlier analysis where values 
of the critical exponents of an Ising system were deter- 
mined from the expansion of the entropy Sjsr(E,M) in 
the vicinity of the critical point of the finite system, situ- 
ated at E = i?jv ;C and M — 0. Here it is our aim to show 
that the non classical infinite lattice exponents can be 
obtained directly from the density of states if the expan- 
sion of the finite system entropy Sn(E, M) is carried out 
at the critical point E = E c<00 and M = of the infinite 
system. (The critical point of the finite or infinite system 
is defined as the value of the energy E — E c ^, where the 
second derivative [8 2 Sn /0M 2 ]e of the entropy taken at 
M = changes its sign.) 

For the purpose of determining the density of states 
Q N (E,M) = exp Sn(E, M) (in units where ks — 1) we 
have developed a highly efficient algorithm which is based 



upon the Broad Histogram mcthodBi. It allows a pre- 
cise and speedy determination of the density of states as 
a function of one, two or even more parameters should 
it be necessary. In the course of these calculations it 
will be demonstrated that very good values for the order 
parameter critical exponent can be obtained by a mi- 
crocanonical analysis already for relatively small system 
sizes. 

The paper is organized in the following way. In the 
next Section we present the numerical method which 
permits us to compute the density of states of classi- 
cal spin systems with very high accuracy. Section III is 
devoted to the analysis of the microcanonically defined 
spontaneous magnetization and to the determination of 
the corresponding critical exponent. A short summary 
concludes the paper. 



II. NUMERICAL METHOD 

One of the aims of computer simulations in statistical 
mechanics is the determination of averages of observables 
X of the many particle system under study: 



(1) 



The space T M contains an enormous number A M 



of 



crostates and the sum in (Q) runs over all of them. W(fi) 
is a probability distribution on that space and X(/i) is 
the value of X in the microstate /i. 

In all cases of interest W(/i) and X(fx) depend on fx 
only via a few parameters P = (P a , Pb, ■■•), i.e., 



and 



W(fi) = W(P{n)) 



X(fi) = X(P(p)) 



(2) 



(3) 



Examples for P_ might be the interaction energy E of a 
magnetic system and its magnetization M. An example 
for W(P) could be the canonical distribution 



W(E, M) 



exp [-/3(E-hM)] 
E exp[-0(E'-hM')Y 

E'M' 



(4) 



The parameters P characterize the macrostate of the sys- 
tem and examples for X(P_) might be the parameters 
themselves, their fluctuations as e.g. (E— < E >w) 2 or 
any other function of them. 
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The unwieldy sum in ([[J) is greatly reduced if the de- 
generacies fi(P) of the macrostates which are character- 
ized by the parameters P are known. Then: 



/ w 



(x) w = Y,n(P)w(P)x(p). 



(5) 



f2(P) is also called the density of states DOS. The number 
of terms in the sum (|l|) is already gigantesque for a system 
of only a few hundred particles: A p is of order f N where 
A" is the number of particles in the system. / is equal 
to 2 for an Ising system and for a g-states Potts model 

/ = «■ _ 

Compared to that, the sum in (0) can be performed 

with very little effort indeed. Its number of terms is only 

of order N v where v is the number of parameters. Once 

f2(P) is known, the sum in (|^) can be performed in the 

twinkle of an eye by any desktop computer. 

Now the difficulty of calculating averages of X has 
been shifted to the determination of the DOS. This does 
not seem to be any easier, but we have developed a 
new algorithm to construct fi(P). Before the algorithm 
is explained in detail we want to point at yet another 
advantage which is offered at no extra cost when the 
density of states Sl(E, M) has been determined. As 
lnfl(E,M) = S{E,M) (setting k B = 1) is the entropy 
as a function of its natural variables, the thermodynamic 
relations can be derived directly from S(E, M) by dif- 
ferentiation. One thereby obtains the microcanonically 
defined temperature, specific heat, susceptibilities etc. 

This is an extra benefit on top of the possibility of hav- 
ing the mean values of any function X(P_) (e.g. X(E, M)) 
and for any desired probability distribution W(P) at 
hand. 

The algorithm is a variant of the transition variable 
methodH~H Its starting point is an extremely simple ob- 
servation: Consider any mechanism which generates a 
new microstate fi when it is applied on the microstate 
fi. An example for such a mechanism could be a sin- 
gle spin Hip in an Ising system. When the single spin 
flip mechanism is applied to all the A" spins of all the 
f2(Pj) microstates which belong to the macrostate which 
is characterized by P_ it then A" f2(Pj) new microstates are 
created. Let Vij be the number of them which belong to 
the macrostate characterized by Pj . As any spin flip can 
be reversed: Vji = Vij . In other words there are as many 
connections (by single spin flip) from level P 4 to level P • 
as there are connections from level Pj to P_ { . Thus if one 
chooses one microstate of level Pj at random, selects one 
spin at random, and flips it, the probability of arriving at 
level Pj will be Vij/(N fi(Pj)) and vice versa if we choose 
one of the spins in one of the states of level Pj at random 
and apply the mechanism to it, then the probability of 
arriving at level Pj is Vij /(N Q(Pj)). 

This simple observation is used to create an algorithm 
which serves two purposes: 

i) It determines the rate of attempts = 
Vij/(ATn(P i )) to go from a starting level Pj to another 



level P_j which can be reached by the chosen mechanism. 

ii) It leads to an equal probability of visiting the levels 
irrespective of their degeneracy. 

Once the Tij have been found it is an easy matter 
to construct the whole surface f2(P) from the ratios 

Ta/T ji = ii(P d )/a(p 4 ). 

In order to calculate the Tjj one writes down the total 
number Zi of applications of the mechanism while in level 
Pj and also the number of applications Z^ which would 
lead to level P- and this in both cases, if the step is 
executed or not. Then Tjj = Z^jZ^. 

As a reversible transition mechanism we choose single 
spin Hips. Before the spin is flipped the system is charac- 
terized by the parameter values Pj. A spin is chosen at 
random. Its Hip would lead to a new state characterized 
by Pj • We add 1 to Z.- L and also to Z^ whether the flip is 
accepted or not. After a large number of attemped spin 
flips ZijjZi approches T^. The demand ii) is fulfilled, 
if steps from i to j and from j to i are executed with 
the same frequency. This can be achieved by introducing 
probabilities for the acceptance of executing a step. This 
probability of accepting the step from i to j is equal to 
1 if < Tji and it is Tji/T^j otherwise. This equalizes 
the probabilities 



p(P.,P l + AP)=p(P l + AP,P t ) 



(0) 



of going from a level Pj to the level Pj + A P and of the 
reverse transition. 

Choosing nonzero, but otherwise arbitrary values for 
Zi and Z^ for all i and j and starting from any one of the 
states of the system one rapidly builds up good estimates 
for the transition rates Tjj. When this point is reached 
all levels Pj are visited with the same frequency and the 
quality of the transition variables is improved with the 
same rate over the whole parameter range. 

The procedure-, [adopted here differs from other 
implementationsffliiy of the transition variable method 
in several respects: The transition variables are updated 
at every single spin flip, in return we do not register all 
the possible transitions which are possible from a given 
microstate, but only those ones which are attempted by 
the algorithm, whether they are successful or not. Fur- 
thermore the transition probabilities Tjj are updated dur- 
ing the whole run and the acceptance rates are governed 
by the most recent value of Tjj from the beginning of the 
simulation to its very end. 

This method gives us the freedom of restricting the 
calculations to a chosen number of channels in energy 
or magnetization direction. Whereas for not too large 
systems all values of the energy and of the magnetiza- 
tion may be allowed in a single run, for larger system 
sizes narrow bands in energy direction covering all mag- 
netization channels may be used. In the limit of only 
one energy channel per bancLoux. approach differs from 
the well-known Q2R methodErtj, as for a fixed energy 
value all the possible magnetizations are visited with the 
same frequency. 
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Our algorithm has been subjected to several stringent 
tests. We have calculated f}jv(P) for a 32 x 32 Ising 
model and we have compared the resulting entropy with 
the exact data of BealeE-X Sn(E) as well as its first 
and second derivate match the exact values almost per- 
fectly. For example, 7 x 10 5 Monte Carlo sweeps yield 
for the inverse microcanonical temperature S' N (E) an av- 
erage error as small as 0.025% when compared to the 
first derivative of the Beale data. In a plot of the spe- 
cific heat c(E) = -(S' N (E)) 2 / S'^^E) versus temperature 
T(E), where 1/T(E) = (3*(E) = S' N (E), one cannot dis- 
tinguish the simultated data from the exact result. For 
the same system we have also calculated Clpf(E, M) from 
which we again find Qn(E), now by summation over M. 
In addition we find Qn(M) by summation over E and 
there is exellent agreement with the exact result which is 
Nl/[ N + M \ N ~ M \]. 

The algorithm compares very favorably to other meth- 
ods proposed for the computation of the density of states. 
Its efficiency for the 32 x 32 Ising model is slightly better 
than that of the multiple-range random walk algorithm 
of Wang and LandauL^I where an average error of 0.035% 
was obtained for 7 x 10 5 Monte Carlo updates. It must 
be stressed that our method is completely different to the 
Wang/Landau method where a rough estimator of the 
DOS is built up rapidly due to a multiplicative update of 
Q(E) every time the energy level E is visited. The error 
introduced by this multiplicative update has to be oliipir 
nated before obtaining the DOS with high precisiont3li3. 
In our present approach, no such error is introduced dur- 
ing the course of the simulation. 

As we will show in the following, the new method en- 
ables us to compute the density of states with a high 
accuracy even for large Ising systems. Furthermore, pre- 
liminary results for complex systems with a rough energy 
landscape are very promisingEj. 

As a first application, we compute in the next Sec- 
tion the microcanonical order parameter in the two- 
and three-dimensional Ising models. It must be noted 
that data of extremely high accuracy, rarely needed in a 
canonical study, are mandatory for a reliable determina- 
tion of microcanonically defined quantities. 

III. MICROCANONICALLY DEFINED 
SPONTANEOUS MAGNETIZATION 

In a microcanonical analysis of finite systems with 
N = L d spins (d being the number of space dimensions) , 
the object of interest is the density of states ilw(E, M) as 
a function of the energy E and the magnetization M. We 
define the spontaneous magnetization M sPt N{E) as the 
value of M where the entropy S N (E, M) = log fl N (P, M) 
at a fixed value of the energy has its maximum with 
respect to M. At energies lower than a critical value 
P c .jv the entropy exhibits two maxima at M = ±M SP! tv- 
On approach to this point E c ^ the order parameter 



M sPt ]\r vanishes with a square root behaviour in the fi- 
nite systemcl. 

In order to achieve very good statistics we have re- 
stricted the calculations to a narrow band which covers 
only five channels in energy direction, but stretches over 
all possible values of the magnetization. All single spin 
flips within this band are allowed. Spin flips which would 
end up outside the band are rejected. We build up the 
Tij and from them construct the DOS in magnetization 
direction. A single flip leads from a microstate in the 
level P_, L — (Ei,Mi) with entropy Stj(Ei,Mi) to a state 
of a neighbouring level P ■ = (Ej , Mj) with Ej = Ei + irj 
and Mj — Mi + 27. Here 7 = ±1, whereas 77 = 0, ±1, 
±2 for a two-dimensional square lattice and r\ = 0, ±1, 
±2, ±3 for a three-dimensional simple cubic lattice. The 
entropy at the macro state at P_j can then be expressed 
as: 

S N {Ej , M 3 ) = S N [Ei , Mi) + 0* {Ei , Mi) 4 r? 

+n(E u M i )2 1 +... (7) 

where /3*(E,M) = 8S N (E,M)/dE is the inverse tem- 
perature and n(E,M) = dS N (E,M)/dM. Because of 
the densely spaced energy and magnetization levels (on 
the intensive scale) higher order terms in (Q) may savely 
be neglected. From ratio of one transition observable 
pair Tij/Tji = n(P,-)M£i) = exp [4?7/3* + 2 7 /j] follows 
one linear combination of (3* and /i. In two dimensions, 
the ten states at the magnetizations M and M + 2 and 
energies E + Aij within the band of states to which the 
calculation was restricted, are connected by 12 transition 
variable pairs. Averaging over the 12 equations which 
determine (3* and /i considerably increases the statistical 
relevance of the data. In three dimensions, the number 
of transition variable pairs is even larger. 

Thus, the transition variable method directly provides 
the derivative /^(P, M) of the entropy which we need for 
the determination of its maximum. 

From the zeros of /x(P, M) the spontaneous magne- 
tization is determined with high precision as shown in 
Figure 1. The variation of the spontaneous magnetiza- 
tion per spin, m sPt N = M SPi n/N, as a function of the 
specific energy, e = E/N, is displayed in Figure 2 for two- 
dimensional Ising models on an infinite and on a small 
square lattice with 32 2 spins. The curve for the infinite 
lattice has-been drawn from the exact results for M(T) 
and P(T)li3. Both curves coincide at low energies. In 
our analysis, M sp {E) for energies lower than the critical 
value P C;00 of the infinite system is fitted by 

M sp {E) = A\e*f s " ff{B) (8) 

with an energy dependent effective exponent /3 £ie //(P). 
Here, e* — (P Ci00 - E)/{E C<00 — E g ) is the reduced en- 
ergy. Eg denotes the ground state energy. The effective 
exponent, which is given by (3 £ ,eff = dlnAf sp /dln |e*|, 
yields the criticaLexponent (3 e in the limit P — ► P c .oo. 
Contrary to Ref.l3 the critical expansion refers to P Ci00 
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and not to the critical energy of the finite system. The 
procedure has been performed for a series of system sizes 
N% < N2 < N3 . . . In the energy range where Ni and 
Ni + i yield the same exponent [3 £ _ e ff(E) this exponent is 
selected for the plot in Figure 3. This value of (3 £je ff(E) 
is used up to the energy where the two exponents start 
to differ from each other, i.e. until finite size effects show 
up. Then we plot the common exponent for the system 
sizes Ni + i and A^ + 2 until these begin to disagree and so 
on. The largest systems for which we have performed 
these calculations were 700 x 700 in two dimensions and 
80 x 80 x 80 in three (these are presumably the largest 
Ising systems for which microcanonical quantities have 
ever been determined) . This approach has the advantage 
that the finite-size effects arc monitored closely. Further- 
more, one avoids to simulate unnecessarily large systems 
at low energy where small systems already yield data es- 
sentially free of finite-size effects. 
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FIG. 1. Derivative fi{E ,M) = 8S(E ,M)/dM of the en- 
tropy of a three-dimensional Ising model with 20 x 20 x 20 
spins computed at the energy Eo/N — —1.2. The sponta- 
neous magnetization M 3P: n/N — ±0.441 is obtained from 
n(Eo,M) = 0. For energies E < E C: n, a local minimum of 
the entropy is observed at M — 0. 

The critical exponent (3 £ is related to the commonly 
used exponent (3 which describes the behavior of the 
order parameter with respect to the temperature by 
/3 £ = f3/(l — a) where a is the usual specific heat crit- 
ical exponentcl. In two dimensions where a = we ex- 
pect j3 £ = 1/8, in three dimensions with (3 = 0.327 and 
a = 0.11 we expect (3 £ = 0.367. 

Also shown in Figure 3 is the corresponding effective 
exponent derived from the exact result M(E) of the in- 
finite two-dimensional Ising modclix Obviously, the nu- 
merical data follow closely the exact result, thus demon- 
strating the reliabilty of our algorithm for computing ac- 
curate microcanonical quantities even for large systems. 




FIG. 2. Microcanonically defined order parameter vs spe- 
cific energy for two-dimensional Ising models defined on an 
infinite square lattice (full line) and on a finite square lattice 
containing 32 x 32 spins (dashed line). The finite-size spon- 
taneous magnetization vanishes at a well-defined finite-size 
critical point. 




FIG. 3. Effective exponent /3 e , e // as function of the reduced 
energy |e*| obtained for the Ising model in two (circles) and 
three (squares) dimensions. Only error bars larger than the 
symbol sizes are shown. The dashed lines and the arrows 
indicate the expected values of the order parameter critical 
exponent j3 e in two and three dimensions. The full line is 
the effective exponent derived from the exact result M(E). 
Note its rapid decrease towards the value 1/8 very close to 
the critical point. 

Figure 3 exhibits two pecularities of the 2D system: 
(i) The strong N dependence of E c ^ makes it extremely 
difficult to obtain data points for |e*| < 0.2. (ii) The 
effective critical exponent a e f / (t) of the specific heat c(t) 
of the infinite system is finite for finite t = (T c — T)/T c , 
but precipitates to zero as t — is approached. This 
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causes the dramatic behavior of (3 E . e ff(e*) = /3 e ///(l — 
a e ff) in the vicinity of |e*| = 0. 

Remarkably, good values of the exponents are already 
found for the smallest systems considered here, namely 
L = 50 in two and L — 10 in three dimensions. This is 
especially true for the three-dimensional case where, over 
the whole energy range considered, the value of the effec- 
tive exponent does not noticeably differ from the value 
of the critical exponent. Even for energies more than 30 
% below the critical energy and for systems consisting 
of only 1000 spins, the obtained effective exponent is al- 
most identical to the critical exponent. This is in marked 
contrast to a similar canonical analysis of the behaviour 
of the order parameter with respect to the temperature 
where large corrections to scaling lead to an effective tem- 
perature dependent exponent which differs by almost 50 
% from the value of the .critical exponent at tempera- 
tures some 30 % below T c iB. Whereas it is obvious from 
Figure 3 that corrections are present when analysing the 
two-dimensional Ising model, these corrections are again 
small compared to the corrections observed for the tem- 
perature dependent effective exponent derived from the 
order parameter M(T). 

The origin of the different behaviour of /3 E , e // as func- 
tion of the reduced energy in two and three dimen- 
sions can be traced back to the logarithmic divergence 
of the specific heat in the two-dimensional Ising model. 
We expect a behaviour similar to that observed in the 
three-dimensional Ising model for all models with a non- 
vanishing specific heat critical exponent. Future investi- 
gations of various models will be needed to clarify this 
point. 
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IV. SUMMARY 



We have developed an algorithm which rapidly builds 
up transition variables for any reversible rule of jump- 
ing from one macrostate of a system, characterized by a 
set P of parameters, to another one. The most recent 
values of the transition variables govern the further path 
of the system through phase space in such a way that 
equally good statistics is obtained over the whole area 
which is covered by the calculation. As an example we 
have chosen the two- and three-dimensional Ising mod- 
els with P = (E, M) together with single spin flips as 
the transition mechanism. From the transition observ- 
ables we obtain highly accurate data for dS/dM which 
allow a precise determination of the microcanonically de- 
fined spontaneous magnetization M sp {E). An analysis of 
M sp (E) for different system sizes yields excellent values 
for the order parameter critical exponent /3 £ . 
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